このチュートリアルでは3次元の画像解析を行ないます。基本的な考え方は2Dの二値画像解析と同じなので、まずそちらのチュートリアルを終わらせてからこちらにトライしてください。
ここでは50x50x50のサイズの画像を解析します。3次元画像の表現方法は色々ありますが、ここでは各層ごとの断面画像を用います。
つまり50x50の画像を50個用意します。
data/ ディレクトリにファイルが置かれています。ファイル名は断面の順に付けておきます。
ls data/
まずはこの 50x50x50 の画像を読み込んで 50x50x50 の白黒配列を作ります。
# 必要なライブラリの読み込み
import imageio
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
3次元配列を作る前に 0000.png を読み込んでどんな画像か確認しましょう。
image_0000 = imageio.imread("data/0000.png")
image_0000.shape, image_0000.dtype, np.max(image_0000), np.min(image_0000)
これは50x50の配列で、0 から255の値を持っています。つまりこの画像はグレイスケール画像のようです。中身を見てから実際に表示してみましょう。
image_0000
plt.imshow(image_0000, "gray")
では50個の画像を読み込んで積み重ねます。
pict = np.stack([
imageio.imread("data/{:04d}.png".format(n)) > 128
for n in range(50)
], axis=0)
pict.shape, pict.dtype
50x50x50の画像が完成しました。ではこれを3次元可視化しましょう。plotly.graph_objectsと
homcloud.plotly_3d モジュールを使います。
詳しくはAPIドキュメントを見てください。
Bitmap3d関数をここでは使います。コメントアウトの所を有効化すると不透明度を変えることができます.
# 必要なモジュールを読み込む
import homcloud.plotly_3d as p3d
import plotly.graph_objects as go
fig = go.Figure(data=[p3d.Bitmap3d(pict)], layout=dict(scene=p3d.SimpleScene()))
# fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.show()
パーシステント図を計算します。hc.PDList.from_bitmap_distance_functionを使います。
2Dの場合と同じです。
import homcloud.interface as hc
pdlist = hc.PDList.from_bitmap_levelset(hc.distance_transform(pict, signed=True), save_to="binary-3d.pdgm")
以下のようにしてファイルを読み込むこともできます。
pdlist = hc.PDList("binary-3d.pdgm")
1次のパーシステント図を調べていきましょう。
pd1 = pdlist.dth_diagram(1)
pd1.histogram((-15.5, 10.5), 26).plot(colorbar={"type": "log"})
(-4, 4) の所に2つbirth-death pairがあるようです。これを確認しましょう
pairs = pd1.pairs_in_rectangle(-4, -4, 4, 4)
pairs
確かに2つあります。
この2つのペアがどのような穴を表現しているのか見てみましょう。1次のパーシステントホモロジーを見ているわけですから、 何らかのリング構造を見ているはずです。それを計算しましょう。
ここでは optimal cycle というものを使います。ポイントクラウドの解析ではoptimal volumeやvolume-optimal cycleというものを用いましたが、
それと似たものです(ちょっと違います)。Pair.optimal_1_cycleで計算します。
optimal_1_cycles = [pair.optimal_1_cycle() for pair in pairs]
まずはこれを可視化しましょう。元のボクセルデータを重ねあわせて表示します。
fig = go.Figure(data=[
cycle.to_plotly3d() for cycle in optimal_1_cycles
] + [
p3d.Bitmap3d(pict, color="black", name="bitmap"),
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.1, selector=dict(name="bitmap"))
fig.show()
回転させたりして観察してください.opacityを変更するなどしてみてみましょう. ボクセルデータの場合には1次のbirth-death pairに対応するのはトンネル構造である、というのが この可視化を見るとわかると思います。
次にこのリング構造がどういうパスを通っているのかを見ます。Optimal1CycleForBitmap.pathを使います。通過点のピクセルの座標が得られます。
optimal_1_cycles[0].path()
optimal_1_cycles[1].path()
pd0 = pdlist.dth_diagram(0)
pd0.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
(-16.0, -10.0), (-15.0, -10.0) という2のbirth-death pairを詳しく調べていきましょう。
ボクセルデータの0次と2次の逆解析にはBitmapPHTrees というクラスを使います(1次には使えないので、上で説明した optimal_1_cycle を代わりに使いました)。
BitmapPHTrees.for_bitmap_distance_functionはPDList.from_bitmap_distance_functionと似たインターフェースで
逆解析用の情報を計算します。
この PHTrees という名前は、ここで計算されるものが木構造であることを意味しています。
n次元のビットマップデータに対して0次とn-1次の2つの情報を同時に計算して、2つの木構造を持っています。
この木構造を取りだすためには PDList.bitmap_phtrees というメソッドを使います。
hc.BitmapPHTrees.for_bitmap_levelset(hc.distance_transform(pict, signed=True),
save_to="binary-3d-tree.pdgm")
以下のようにして計算したデータを読み込みます。
pdlist_with_tree = hc.PDList("binary-3d-tree.pdgm")
0次の情報を取り出すのは bitmap_phtrees(0) を使います。
phtree_0 = pdlist_with_tree.bitmap_phtrees(0)
(-16.0, -10.0), (-15.0, -10.0)という2つのbirth-death pairの情報を取り出すのには
pair_nodes_in_rectangleを使いましょう。これは指定した範囲内にある birth-death pair (に対応する木構造のnode)をすべて取り出します。
nodes_0 = phtree_0.pair_nodes_in_rectangle(-16, -15, -10, -10)
nodes_0
取り出したデータは BitmapPHTrees.Node というクラスのインスタンスです。次にこれらの情報を可視化しましょう。
fig = go.Figure(data=[
n.to_plotly3d() for n in nodes_0
] + [
p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()
0次のPHなので,ボクセルの内側の領域に 対応しています.実はこの2つのbirth-death pairは木構造上の親子関係にあり,一方が他方の部分集合になっています.よく見ると赤い領域の一部が青くなっているのが見えると思います.
pd2= pdlist.dth_diagram(2)
pd2.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
(4, 6)というbirth-death pairを調べます。逆解析には
bitmap_ph_trees(2)を使います。2が次数です.
phtree_2 = pdlist_with_tree.bitmap_phtrees(2)
(4, 6) という birth-death pairを調べましょう。まず pair_nodes_in_rectangle で (4, 6) という birth-death pair (に対応した tree の node) をすべて取り出します。
nodes_2 = phtree_2.pair_nodes_in_rectangle(4, 4, 6, 6)
nodes_2
2個の 以下が可視化した結果です。くりぬかれた領域の奥のほうがこのbirth-death pairに対応していることがわかると思います。
fig = go.Figure(data=[
n.to_plotly3d() for n in nodes_2
] + [
p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()
volumeメソッドを使うことでこの領域の各座標を得ることもできます。
nodes_2[0].volume()